Very large scale simulations of the RSOS model in four dimensions 
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We study the Restricted Solid on Solid (RSOS) model for surface growth in spatial dimension d = 4 by 
means of a multi-surface coding technique that allows to analyze samples of unprecedented size. For such 
large systems we are able to achieve a controlled asymptotic regime where the typical scale of the fluctuations 
are larger than the lattice spacing used in the simulations. A careful finite-size scaling analysis of the critical 
exponents clearly indicate that d = 4 is not the upper critical dimension of the model. 
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The Kardar-Parisi-Zhang (KPZ) equation fTl is possibly the simplest, most studied, and yet not fully understood model for 
out-of-equilibrium surface growth. The equation describes the time evolution of the height h{r, t) of an interface above a 
d— dimensional substrate; 

dth{r,t) = :^V^h{r,t) + ^\Vh{rA)\^ + T^{r,t) , (1) 

where is the diffusion coefficient, A is the strength of the non-linear growth rate term which is responsible for the h ^ ~h 
symmetry breaking with respect to the growing direction, and r/(r, t) is a Gaussian white noise of amplitude D: 

(r/) = , (?7(r, t)77(r', t')> = 2i?<5'^(r - r')S{t - t') . (2) 

The KPZ equation describe many relevant growth processes, such as the Eden model, ballistic deposition, interface growth in 
disordered medium. It is also related to many other physical phenomena apparently unrelated to surface growth, such as Burgers 
turbulence, dynamics directed polymers in random media, dissipative transport in the driven-diffusion equation |2|. 

The scaling properties of the height's fluctuations W2{L, t) — (/i^(r, t))^ — {h{r, t))1 (with the notation (. . . )r we indicate a 
spatial average over a macroscopic hypercubic box of linear size L over the d— dimensional substrate) characterize the univer- 
sality class of the model. More precisely, for a system of size L, W2{L,t) ^ L'^^f{t/L'^), where the scaling function is such 
that f{x) — > const for x oo and f{x) ^ x^^l^ for a; — ?> 0. The peculiar behavior of / imply that W2{L, t) ^ L"^^ for t^L^ 
and W2{L, t) ^ t^^l'^ for t L^. Due to an infinitesimal tilt symmetry of Eq. >/i + r-e, r— >-r — Ate), the two critical 

exponents are related by the scaling relation x + ^ = 2, which is believed to be valid at any dimension d f2\. 

A complete understanding of Eq. ([TJ, and in particular the determination of the two critical exponents x, z for any spatial 
dimension d > 1 (at d = la fluctuation-dissipation theorem leads to the exact result % = 1/2, z = 3/2) , turns out to be 
extremely difficult for two main reasons: (i) we are dealing with an intrinsically out-of-equilibrium phenomenon where the 
standard equilibrium toolbox must be used with care, (ii) perturbative renormalization schemes are not adequate for describing 
the strong coupling regime (i.e. where the parameter A is relevant). 

The existence of an upper critical dimension d^, i.e. the substrate dimensionality d above which the fluctuation of the model 
become irrelevant (x = 0), is one of the most controversial unsolved theoretical issues related with Eq. ([T}. The determination 
of d„ would be a most relevant achievement since, as customary in equilibrium critical phenomena, its knowledge constitutes 
the first step for a controlled perturbative expansion around it. The quest for du has been around for more than twenty years 
OmSil and the different predictions range from du ~ 2.8 to d^ = oo. Analytical estimates using the mode-coupling theory yield 
exact results in rf = 1 1T61 . Their extension to higher dimensions hints for a c?u = 4 under different self-consistency schemes 
lISJH. The same value for du is also supported by different field-theoretic approaches ll3l [T0l - [T2l . 

At odd with what predicted by the previously mentioned field-theoretic approaches, both direct numerical integration of KPZ 
equation IflTl . and simulation of systems belonging to the KPZ universality class lfT3l[T4l[T8l[T9l indicate that du > 4, while the 
real-space renormalization group approach 1 15| predicts du — oo. 

Such a long standing controversy is the consequence of the difficulties inherent to both analytical and numerical approaches. 
Most of the assumptions made on the functional structure of the sought solution in the different field-theoretic analysis, as well 
as the approximations made in the mode-coupling theories are, in general, not completely under control. On the numerical side 
the most severe problem is due to the fact that simulations in high spatial dimensions d > 4 are computationally very heavy, and 
the systems under analysis must be limited in size. As a consequence, the different fitting procedures must deal with controlled 
finite-size scaling procedures to yield reliable estimates of the critical exponents. Under this perspective, particularly relevant 
is the observation that for lattice models in the KPZ universality class, a controlled asymptotic regime is achieved only when 
typical scale of the fluctuations is larger than the lattice spacing used in the simulations or, more precisely, for u'2 > 1 |9|. The 
former inequality is very stringent from the computational point of view since it requires very large lattices to be fulfilled: the 
estimates presented in 1 14] suggest indeed that for the 4— dimensional RSOS model, to which this letter is addressed, the W2 > 1 
inequality starts being verified for lattice size larger then L w 32, whereas the larger system size analyzed in the same paper was 
L — 28 that, at the best of our knowledge, remains the larger system in 4 dimensions simulated so far. 

To settle the controversy, at least in the 4— dimensional RSOS model case, we decided to analyze samples of unprecedented 
size: we have been able investigate the steady state scaling regime t ^ for lattice size volumes up to = 128^ = 268435456 
sites, a factor 437 off the largest simulation known in literature 1141 . We can also study the dynamic scaling regime of lattices 
of V — 256* — 4294967296 sites, but for such a large size we have been able to investigate the asymptotic scaling regime for 
just three samples, due to limitations in our computing facility (most of the data for L = 256 are at not too large t and they have 
been used only in fig. (1): they appear only in the region t/L^ < 8). 

The RSOS can be simulated in the following way: at any time t we randomly select a site i on the d— dimensional lattice and 
we let the surface height hi at that point to grow of a unit hi{t + 1) = hi{t) + 1 only if maxj^Qi \hi{t) — hj{t)\ < 1, being di 
the set of 8 nearest neighbors of i in d = 4 (note that we will assume periodic boundary conditions). 
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We simulated RSOS growth using two different algorithms based on a very efficient multi-spin coding technique ll20ll origi- 
nally developed for disordered spin system, and later generalized to deal with the RSOS model 1141 : 

• Multi Surface (MS) coding: we can simulate, with basically the same cost of one single surface simulation, Ni, copies 
of the system, Nt, being the number of bits in an computer world (usually 32, 64, 128 and 256). We transform the basic 
operations (like summing spins for computing the effective force) into Boolean operations, and we exploite the fact that 
when, for instance, the computer is calculating an AND logical bit, it is indeed doing that operation Ni, times at once, i.e. 
for all the bits of the world. Unfortunately computational efficacy of this method is counterbalanced by the memory load, 
making it unpractical for analyzing samples of linear size larger than L — 64, at least on the computers we have access to. 

• Multi Lattice-site (ML) coding: for sample of linear size L = 128, 256 we have developed a new multi-spin coding 
representation in which a single surface at time is simulated, but we lump together Ni, height-difference local variable in 
a single computer word of Ni, bits. We will refer to this second method as the Multiple Lattice-site (ML) algorithm (see 
tablejl]). In this algorithm in the first half step we update the even spins (i.e. the spins where ix + iy + iz + H is even) 
applying the standard RSOS procedure with probability 1/2 to each spin (with probability 1/2 the spin is not updated), in 
the second half step we apply the algorithm to odd spins. Some programming care must be used with periodic boundary 
conditions: in the simplest version we have used, L must be an even multiple of Ni,. Moreover it is crucial to use a good 
random number generator, where all the bits are random. 



L 


sweeps 


samples 


type 


time (h.) 


8 


524000 


1024 


MS 


4 


16 


524000 


1024 


MS 


6 


31 


524000 


1024 


MS 


121 


32 


524000 


1024 


MS 


139 


33 


524000 


1024 


MS 


158 


64 


131000 


512 


MS 


5376 


128 


512000 


64 


ML 


7680 


256 


130000 


64 


ML 


504 



TABLE I. In this table we display the lattice linear size L, the number of montecarlo sweeps (full lattice update), the number of samples and 
the simulation type (MS = multi-surface coding, ML = multi-lattice coding), and the overall computational time in hours. 



We simulate 4— dimensional lattices of volume V — L"^ for lattices of linear size L = 8, 16, 31, 32, 33, 64, 128. For the two 
largest lattice (L = 64, 128) we run the ML algorithm, while for the rest we run the MS algorithm. We decided to consider 
L = 31, 33 for checking that there are no periodicity issues with the random number generator A summary of our simulations 
is provided in table |l] 





X 




A2 


B2 


Az 




Ai 


54 


NEW 


0.2537(8) 


1.11(9) 


0.171(1) 


0.37(6) 


0.0319(3) 


-1.0(2) 


0.100(1) 


0.38(8) 


OLD 


0.255(3) 


0.98(9) 


0.170(1) 


0.37(3) 


0.0321(2) 


-0.7(1) 


0.100(1) 


0.46(4) 



TABLE II. In this table we display the best fit values together with their statistical error of the parameters defined in Eq. |4]). The first row 
refers to the actual data presented in this work, the second is taken from fT4l 



The simulations aim at achieving a fair sampling of the asymptotic regime. To do so, at any time t and for each sample (both 
in ML and MS type of simulation) we evaluate the first three connected moments Wn{L, t) — ^ (^(^)))"/^' where 

{h{t)) = X^iLi hi{t) /V, and n — 2,3, 4. We thus define our asymptotic (in time) estimate as: 

1 

Wn{L)^- rj. , , Y] Wn{L,t) . (3) 

J - 1 + i 

t=l 1 

We are careful to choose both To, and Ti — To large enough to guarantee: (i) convergence to the asymptotic regime, i.e. that 
Ti ^ L^, (ii) a large enough sampling of statistically uncorrelated measures of w„(T, t). In practice we consider consecutive 
measuring windows of length 64, 128, Ti, Tq, so that Tq is the last measure in the simulation and Ti = To/2. This choice is 
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very conservative, since eventually we use just the second half of the simulation, but at the same time it allows us to check with 
very high reliability whether or not we have reached the asymptotic state: a quick look at Fig. [T] will comfort our confidence of 
having chosen, even for the larger sample size, a Tq ^ L^, at least of a factor 100 off. We are now ready to determine the critical 
exponents of the asymptotic behavior of w„, which scales as L'"'^ at the leading order, by fitting simultaneously the following 
moments: 

W2= A2L''^{l + B2L-^) 

w^^A^L^^il + B^L-^) (4) 

where cj is the leading finite size scaling exponent. The fit involves the simultaneous determination of 8 parameters whose 
best-fit value is reported in the first row of table|ll](the fit yields a chi-squared root mean square deviation of 2.05). Interestingly 
enough the fit presented in fT?] agrees very well with the new data as we can clearly appreciate qualitatively in Fig.|2j and more 
quantitatively by comparing the two rows in table [ll] With respect to the ?z;2 > 1 inequality issue, a glance at Fig. [2] shows 
unambiguously that the scale of the typical fluctuations, for all lattice size larger than L = 31, verify the inequality. We do not 
see any change in the scaling behavior of the three cumulants around the cross-over region L w 30, moreover, the fact that the 
old fit presented in [141 (in a regime ■W2 < 1) agrees so well with our larger lattice size simulation (see again Fig.|2| indicate 
that the simulations performed for L < 28 were able to capture fairly the asymptotic scaling regime. To see more clearly the 
finite-size corrections of x we determined the effective exponent Xr!^ as the discretized logarithmic derivative of Eqs. ^ which 
in our case reads: 



xfiL)^ rTr/ (5) 



where L/L' — 2, and n = 1, 2, 3. In Fig.jsjwe display x'^ a function of (note that we discarded for the sake of clarity the 
L = 31, 33 results) for the three cumulants together with the best-fit curves. The fit yields the following results for the critical 
exponents: x ~ 0.2532(5) and uj = 1.14(5) (see also table [ll|. 

A recent has investigated a model of direct polymers in random medium that should belong to the same universality class fTSl. 
In this model one can define an exponent C, that according to the theoretical expectations should be given 



Their results {C, slightly larger than 0.57), is well consistent with our prediction C, = 0.5725(2). 

The numerical technique we have introduced has allowed us to run very precise numerical simulations of the RSOS model in 
d = 4 on unprecedented system size with a limited amount of computational time. Thanks to the accuracy of the simultaneous 
measurement of the three cumulants, the claim that d = 4 is the upper critical dimension for systems in the KPZ universality 
class has to be rejected. Moreover, the typical fluctuation's length-scale of our simulations on samples of linear size L = 128 
and L = 256 are larger than the lattice spacing, and this is a clear indication that: (i) the system reached a controlled scaling 
regime, (ii) the measured scaling exponents are reliable and not affected by a pre-asymptotic cross-over regime. 

There is still a remote possibility that our data are consistent with an upper critical dimension (i„ = 4 of the KPZ equation if 
we drop the hypothesis that RSOS in d = 4 belongs to the KPZ universality class. However, apart from some work in the past 
ETl . this hypothesis does not seem to have support in the mainstream literature on KPZ. 

We thank Enzo Marinari, and Massimo Bemaschi for interesting discussions, Moshe Schwartz for relevant correspondence, 
and Marco Zamparo for reading the manuscript. The numerical simulations presented here were run using facility of the Human 
Genetics Foundation. The European Research Council has provided financial support through ERC grant agreement no. 247328. 
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FIG. 1 . Scaling plot of the rescaled second moment ■W2 / {A2L'^^ {1 + B2L ) ) vs. the rescaled time t/ . 



6 



16 
14 
12 
10 



W2 new data ^- 
W2 old fit 
W3 new data •- 



W4 new data 
W4 old-fit 




60 80 100 120 140 
L 



FIG. 2. The quantities 102 , W3 W4, are displayed as a function of L. Dots with error bars are values obtained by simulations, while solid lines 
are the 8-parameters best-fit reported in 1 14|. The vertical arrow arrow at L = 28 represents the largest size simulated in 1141 . 




FIG. 3. Local slopes of u)2,3,4 are displayed as a function of L ^. Dots with error bars are values obtained by simulations, while lines are the 
8-parameters best-fit reported in table|ll] The solid horizontal line is at x = 0.2538, i.e. the best-fit prediction for the wandering exponent. 



